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Abstract 

In Evans function computations of the spectra of asymptotically constant-coefficient 
linear operators, a basic issue is the efficient and numerically stable computation 
of subspaces evolving according to the associated eigenvalue ODE. For small sys- 
tems, a fast, shooting algorithm may be obtained by representing subspaces as sin- 
gle exterior products [2,8,9,10,6]. For large systems, however, the dimension of the 
exterior-product space quickly becomes prohibitive, growing as (^), where n is the 
dimension of the system written as a first-order ODE and k (typically ~ n/2) is the 
dimension of the subspace. We resolve this difficulty by the introduction of a sim- 
ple polar coordinate algorithm representing "pure" (monomial) products as scalar 
multiples of orthonormal bases, for which the angular equation is a numerically op- 
timized version of the continuous orthogonalization method of Drury-Davey [13,14] 
and the radial equation is evaluable by quadrature. Notably, the polar-coordinate 
method preserves the important property of analyticity with respect to parameters. 
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1 Introduction 



A useful tool in the study of stability of travehng waves is the Evans func- 
tion, an analytic function whose zeroes correspond to the eigenvalues of the 
linearized operator about the wave. More generally, let L be a linear differ- 
ential operator with asymptotically constant coefficients along some preferred 
spatial direction x, and suppose that the eigenvalue equation 

(L-X)w^O (1) 

may be expressed as a first-order ODE in an appropriate phase space: 

= A{x, X)W, hm A{x, A) = A±(A), (2) 

with A analytic in A as a function from C to C^(]R, C"^") and the dimension k 
of the stable subspace 5*+ of A^ and dimension n — A; of the unstable subspace 
U- of A_ summing to the dimension n of the entire phase space. 

Then, one may construct analytic bases of solutions of (1), say wf, . . . 
and Wfc+i, . . . , respectively, spanning the manifolds of solutions decaying 
as a; — > -t-cxo and — cxd by essentially initializing them at infinity with values 
from the stable (resp. unstable) subspace of A_^ (resp. A_) and solving toward 
X — using the ODE (2). The Evans function is then defined as 

- det ...W^ W^.,, ■ ■ ■ W-) , (3) 

where each Wi is the solution of (2) corresponding to Wi] for details, see, 
e.g., [1,20,15,29,30] and references therein. Analogous to the characteristic 
polynomial for a finite-dimensional operator, D{-) is analytic in A with zeroes 
corresponding in both location and multiplicity to the eigenvalues of the linear 
operator L. 

Numerical approximation of the Evans function breaks into two problems: the 
computation of analytic bases for stable (resp. unstable) subspaces of A^ (resp. 
^4-) and the solution of ODE (2) on a sufficiently large interval x e [0, M] 
(resp. x G [— M, 0]). In both steps, it is desirable to preserve the fundamental 
property of analyticity in A, which is extremely useful in computing roots by 
winding number or other topological considerations. 

A difficulty in the second problem is numerical stiffness for k, n — k 7^ 1, due to 
the need to resolve modes of different exponential decay (resp. growth) rates. 
This may be overcome in elegant fashion by working in the exterior product 

space A ■ ■ ■ A G c(fe) (resp. W^_^^ A ■ ■ ■ A G c("-'=)), for which the 
desired subspace appears as a single, maximally stable (resp. unstable) mode. 
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the Evans determinant then being recovered through the isomorphism 



det iw,+ ■■■ Wt W, 



k+l 




W+A---AW+AW,-^,A---AW-. 



The first instance of this "exterior-product method" in the Evans function 
context seems to be a computation carried out by Pego in the Appendix of 
[2]. The method was subsequently rediscovered and further developed by Brin 
et al. [8,9,10] and, independently, by Bridges et al. [6,3]. See also, the earlier 
"compound-matrix method" introduced by Gilbert and Backus [16] and also 
Ng and Reid [24,25,26,27] for the numerical solution of stiff ODE, of which it 
may be regarded as a coordinate-independent implementation. 

The computation of an initializing analytic basis at plus (resp. minus) spatial 
infinity is likewise straightforward in the exterior-product framework, since 
it reduces to the calculation of a single eigenvector. Two quite satisfactory 

approaches to this problem were given in [10] and [6], each of order (?) 




equal to the cost of a matrix inversion or the multiplication of two matrices 
in dimension (j^ x negligible compared to the cost of integrating the 
exterior- product version of (2). 

Together, these two steps give an extremely fast and well-conditioned shoot- 
ing algorithm for the computation of the Evans function, for small values of 
n. However, for equations of large dimension n, such as those that arise in 
complicated physical systems or through transverse discretization of a multi- 
dimensional problem on a cylindrical domain [23] , the exterior- product method 
quickly becomes impractical, since the typical value /c ~ n/2 leads to a working 
dimension ^^^2) growing as n"/^. For example, for the typical values n ~ 10^ 
found in [23], this is clearly out of computational range. The development 
of new numerical methods suitable for stability analysis of large systems was 
cited in a recent A.I.M. workshop on Stability Criteria for Multi-Dimensional 
Waves and Patterns, May 2005, as one of three overarching problems facing 
the traveling wave community in the next generation [19]. 

Discussion of this problem has so far centered mostly on boundary- value meth- 
ods. For example, one may always abandon the Evans function formulation 
and go back to direct discretization/Galerkin techniques, hoping to optimize 
perhaps by multi-pole type expansions on a problem-specific basis. However, 
this ignores the useful structure, and associated dimensionality reduction, en- 
coded by existence of the Evans function. 

Alternatively, Sandstede [28] has suggested to work within the Evans function 
formulation, but, in place of the high-dimensional shooting methods described 
above, to recast (2) as a boundary- value problem with appropriate projective 
boundary conditions, which may be solved in the original space C" for individ- 
ual modes by robust and highly- accurate collocation/continuation techniques. 
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This would reduce the cost to polynomial order C{n)kn'^, where C(-) counts 
the number of mesh points times evaluations per mesh required for system size 
n. He points out, further, that if one is interested only in zeroes of D{-) and 
not component subspaces, then the cost may be reduced by a further factor 
of /c by a root-finding algorithm computing only a single candidate eigenfunc- 
tion in each of the unstable (resp. stable) subspaces (the "bordered matrix" 
method as described in [4,22]). In [18], there were presented correspondingly 
efficient 0{n^) initialization routines prescribing analytic basis/projection for 
the stable (resp. unstable) subspace of (resp. A_), also in the original space 
without reference to exterior products. 

However, up to now, it is not known how to recover analyticity of the Evans 
function in numerically well-conditioned fashion by the above-described collo- 
cation methods. The reason is precisely the (normally advantageous) property 
that errors are uniformly spatially distributed for such methods, whence the 
relative error near spatial infinity, where solutions exponentially decay, is pro- 
hibitively large to track dependence on initializing conditions. Likewise, the 
appealing simplicity/ease of coding of shooting methods is lost in this ap- 
proach. 

Motivated by these circumstances, we introduce in this note an alternative, 
shooting method designed for large systems, couched like collocation meth- 
ods in the original (relatively) lower-dimensional space C" but preserving the 
useful properties of analyticity, simplicity, and good numerical conditioning 
enjoyed by the exterior-product method. Indeed, being based on standard ma- 
trix operations, our method is substantially easier than the exterior-product 
method to code, and has an inherent parallel structure that may be exploited 
in a transparent fashion by the use of a numerical package incorporating paral- 
lel matrix routines; likewise, because it is carried out in minimal coordinates, 
there is no need to take advantage of sparse matrices such as occur in the 
exterior-product method. 

The basis of the new method is to represent the exterior products of the 

columns of W± in "polar coordinates" where the columns of £ 

are orthonormal bases of the subspaces spanned 

by the columns of W+ := {w+ . . . W^+^ and W_ := {w^^^ ■ ■ ■ W") , 

defined as in (3), i.e., W+ = ^+0i+, VF_ = fi_Q;_, and 7± := detQ;±, so that 

where O-^ denotes the jth column of and likewise 

The idea is that the projectivized, angular fiow in Vl should remain numerically 
well-conditioned since orthonormality is enforced, whereas the radial equation. 
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being scalar and linear, is automatically well- conditioned and evaluable by 
simple quadrature. Indeed, this turns out to be the described in the 

remainder of the paper. 

Just as the exterior-product method has ties to the much earlier compound- 
matrix method, our polar-coordinate approach has ties to another well-known 
method, that of "continuous orthogonalization" , introduced by Drury [14] and 
Davey [13] as an alternative to the compound-matrix technique. Specifically, 
our angular equation in principle computes the same Q computed in the con- 
tinuous orthogonalization method. The new feature of our method is the radial 
equation, which restores the important property of analyticity with respect to 
parameters with no loss of numerical conditioning. This represents a subtle 
but important departure in point of view, in that we truly compute exterior 
products and not subspaces or bases thereof, as in past interpretations of the 
(standard) continuous orthogonalization method [5]. As discussed in Section 
3.1, this gives us considerable flexibility in choosing a numerically optimal 
implementation. For the examples considered, we found excellent results (and 
optimal computation cost among methods considered) with the original equa- 
tion of Drury [14]; see Section 4. However, this flexibility might be useful in 
more numerically challenging situations. 

We point out that Bridges [5] has developed a clever "biorthogonal" variant 
of Davey's continuous orthogonalization method that also preserves analyt- 
icity (see Remark 9); indeed, not only the minors of Q but also individual 
columns vary analytically with respect to a parameter. However, this method 
does not preserve orthogonality but only biorthogonality = Ikxk with a 
simultaneously computed "left basis" Q, hence, in addition to requiring twice 
the computation time due to the doubled variable {Q. Q), appears to be inher- 
ently less stable than the other methods considered here. As far as we know, 
this method has never been implemented numerically; it was presented in [5] 
as an interesting "dual" version of the exterior-product method. 

Our numerical experiments indicate that the polar-coordinate method is quite 
competitive in mesh-size requirement with the exterior-product method. Thus, 
the break-even dimension for the polar coordinate vs. exterior product method, 
taking into account competing effects of nonlinear function calls vs. higher 
dimension, seems to be about n > 6 (n > 8 for optimized exterior product 
with sparse matrix solver, which at the moment does not exist). Detailed 
comparisons are given in Section 4. As compared to collocation/continuation 
methods, we expect as for any shooting method that there is a transition 
size above which the stability advantages of collocation outweigh the speed 
and ease of computation of our algorithm. In particular, for "medium-sized" 
systems such as occur in large but genuinely one-dimensional systems such as 
magnetohydrodynamics (n = 15) or combustion with many species (n ~ 10 or 
10^), we expect (though, so far, no study has been made for large systems with 
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either method) that our algorithm is at least competitive with the present 
alternatives for root-finding. And, for the moment, it is the only choice for 
calculating winding numbers in medium-sized systems or above. 

Finally, we point out that, for ultra-large systems for which shooting may 
not be well-advised, our algorithm may equally well serve as the basis of 
analyticity-preserving collocation methods, since uniform errors in orthonor- 
mal bases give good tracking of subspaces along the whole real line. Thus, it 
seems perhaps useful in this context as well. 



2 The algorithm 

2. 1 Derivation 

Our starting point is a comment by Chris Jones [19] that the representation 
in the exterior- product method of products Wi/\ - ■ ■ /\ Wk as the direct sum of 
products of standard basis elements is extremely inefficient, since the less than 
k X n-dimensional manifold of "pure" (i.e., monomial) products is quite small 
in the (^^'^ -dimensional space of direct sums, and one should therefore seek 
more efficient coordinates for the computation. A natural choice is to represent 
exterior products A in "polar coordinates" as (7,$^), where radius 7 e C is a 
complex scalar and angle fl e C"^'^ is matrix of orthonormal column vectors 
Q*Q = /fcxfe whose columns span the subspace spanned by the factors of F, with 
7 chosen so that the product of 7 with the exterior product of the columns 
of fl is equal to A. This representation is unique modulo transformations 
(7, n) {j/ detU,flU), U e C'^^'^ is unitary. The set of angles Q may be 
recognized as the Stiefel manifold described in [4] , a standard coordinatization 
of the Grasmannian manifold of linear /c-dimensional subspaces over C". 

In these coordinates, our computations have a concrete, linear-algebraic in- 
terpretation in which no reference to exterior products appears. Hereafter, 
let 5RM := (1/2) (M + M*) denote the Hermitian part of a matrix M and 

'^M := (1/2)(M — M*) the skew-Hermitian part, and ' denote d/dx. De- 
noting by W+ e C"""^ and W' G the matrices {W^ , . . . ,W^) and 
(PF^i, . . . , FF~) from whose columns the Evans function is determined by 
(3), we have 

W+ = Q^a;+; det a+ = 7+, (4) 
and likewise for W~ . Thus, (3) becomes simply 

D(A) = 7+7- det(Q+, Q-),=o. (5) 
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Fix the Q-evolution by the choice 



= 0, (6) 

removing the ambiguity in representation. This may be achieved by a unitary 
transformation il^m^U e C^^'^ satisfying U' = -n*n'U. (Note that -iim 
is skew-symmetric, since = — 2K(fi*f2a.).) This choice is optimal in 

the sense of minimizing arc length in the Stiefel coordinates, as discussed in 
[4]. It may also be recognized as the one prescribed by a standard continuation 
algorithm of Kato [21] with the orthogonal projection P(A) := flfl*; see [18] 
or Remark 4 below. 

Comparing W = AW = AQa and W = (Qa)' = Q'a + Qa', we obtain 

AQa = n'a + Qa'. (7) 
Multiplying on the left by fl* and invoking (6) and — I, the key equation 

a' = {n*An)a (8) 
relates the two coordinatizations of the desired subspace. 

By Abel's equation, we obtain from (8) and definition 7 deta 

i = trace (Q*AQ)7. (9) 

Finally, substituting (8) into (7), multiplying on the right by a~^, and rear- 
ranging, we obtain 

= (7 - nn*)An, (lo) 

completing our description of the polar-coordinate flow. As described in the 
introduction, ODE (10) is exactly the continuous orthogonalization method 
of Drury [14]. 

Remcirk 1 Equation (9) gives another sense in which is a minimal, or 

"neutral" choice of basis; namely, it is the unique choice for which the gen- 
eralized Abel formula (9) holds up to complex phase. (It holds in modulus for 
any orthonormal basis choice.) 

Rescaled radial flow. Since we ultimately evaluate 7 at x = 0, we may 

strategically introduce, similarly as in [8,9,10] for the exterior-product method, 
the rescaled variables 

7±(x) := 7±(x)e-*'^^^" (n*Aa)±x ^^^^ 

for which the flow near x — ±00 is asymptotically trivial. Our complete algo- 
rithm thus becomes 

n' ^{i -Qn*)A{x,x)fi 

f^trace (n*(A- A_)n)7, ^^^^ 
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with 



L>(A)=7+7-det(0+,0-),=o- 



(13) 



Summarizing, we have 

Proposition 2 For any choice of ^{—oo) and f2(— oo)) with columns span- 
ning the unstable subspace f/+ of A_ := limx^_^A{x,X), there are unique 
solutions of (8), (12) such that 

hm il(x) — Jl(— oo), 

X—^—00 

^lim^7(x) = 7(-oo), (14) 
hm a(x)e-''^'' (n*Af2(-oo))x ^ ^(_oo)/, 

x—*—oo 

and these satisfy W — Q,a, 7 = e"*""^^ {n*Aa,{-oo))x where W is a solution 

ofW'^AW. 

PROOF. Standard asymptotic ODE theory [12] and the above calculations 
relating W, (7, Jl), and 7. 

2.2 Initialization at infinity 

To complete the description of our method, it remains to prescribe the initializ- 
ing values Q(±oo), 7(±oo) in Proposition 2, taking care to preserve analyticity 
with respect to A. Recall the following standard result of Kato. 

Proposition 3 ([21, § II. 4. 2]) Let P{X) be an analytically varying projec- 
tion on a simply connected domain Q. Then, the linear analytic ODE 

r; = [P',Ph; r,(Ao) = rJ (15) 

defines a global analytically varying basis {rj{X)} of the associated invariant 
subspace Range P(A), where denotes d/dX and [A,B] :— AB — BA the 
commutator of A and B. More generally, 

S'=[P',P]S; S{Xo)=I. (16) 

defines a globally analytic coordinate change such that 

S'^PS = constant = Pq- (17) 

PROOF. Relation (17) follows from {S~^PSy = 0, and may be established 
directly by using the key relations PP'P = and (/ - P)(/ - P)'(^ - ^) = 0, 
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which in turn follow by differentiation of the projective equation — P — 0. 
Observing that satisfies the "transpose" ODE, 

(s-^y = -s-^s's-^ = -s-^[p', p]ss-^ = -s-^[p', p], (18) 

we have that both S and its inverse satisfy linear analytic ODEs, hence have 
global bounded analytic solutions in fl by standard analytic ODE theory [12]. 
Finally, Range P = jS" Range Pq is spanned by the columns of SRq, where the 
columns of Rq are chosen to span Range Pq, verifying (15). 

Remcirk 4 Let R{X) — (ri • ■ • rj^ £ C"^'^, be the matrix of basis vectors of 

Range P defined by ODE (15) and L e C''"^" be the matrix whose rows form 
the dual basis o/RangcP*, LR = Jfc G C'''^'^. Then (see Proposition 2.5, [18]), 
the flow (15) is uniquely determined by the property 

LR' = L'R = 0. (19) 

Remark 5 From (18) we see that S (hence R) is unitary if P is self-adjoint 
(i.e., orthogonal projection), since in that case [P',P] = —[P',P], so that 
and S* satisfy the same ODE with same initial conditions I. Likewise, 
the relation P — SPoS~^ — SPqS* shows that S is unitary only if P is self- 
adjoint with respect to some fixed coordinate system ( any coordinates for which 
Pq is self- adjoint). 

Proposition 3 describes a method to generate a globally analytic matrix ^^"(A) 
(the matrix P(A) above) with columns spanning the unstable subspace U~ 
of A- (A) and similarly for S~^, A^. Efficient numerical implementations are 
described in [18]. 

Likewise, we may efficiently compute a matrix Q(— oo. A) at each A whose 
columns form an orthonormal basis for U, for example by either the SVD or 
the QR-decomposition. This need not even be continuous with respect to A. 
Equating 

Q(-oo, A)a(-oo, A) = ly-(A), 

we obtain 

q;(-oo. A) = Q*(-oo, X)W-{X) 

and therefore the exterior product of the columns of W~ is equal to the exterior 
product of the columns of Q{—oo) times 

7(-oo) :=deta(-oo,A) = det(n*(-oo,A)W^-(A)). (20) 

That is, the exterior product represented by polar coordinates (7, Q){—oo, A), 
with 7(— 00, A) defined as in (20), is the same as the exterior product of the 
columns of W~{X) appearing in the definition of the Evans function, in par- 
ticular, analytic with respect to A. 
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With these initiaUzing values (7, Q){—oo, A), we may then efficiently solve (12) 
for Q from x = ±00 to a; = to obtain f2±(0) using any reasonable adaptive 
Runga-Kutta solver, with good numerical conditioning. (The reason, similarly 
as for the exterior-product method, is that n(±oo) is now an attractor for 
the flow in the direction we are integrating; see discussion, [1,15,8,9,10,6].) 
Combining, we obtain D{X) through formula (13). 



3 Further elaborations 

3.1 Numerical stabilization 

We briefly describe various alternatives to the basic scheme (12), designed to 
improve numerical stability in sensitive situations. In the examples we consid- 
ered, such additional stabilization was not necessary. 

Angle equation. It was reported soon after its introduction that the basic 
continuous orthogonalization method (10) can in some situations suffer from 
numerical instability, by Davey [13], who suggested as a variant the Davey (or 
generalized inverse) method 

= (/ - QQ+)AQ, (21) 

where := (i7*i7)~^i7* denotes the generalized inverse. The point is that 
3?(r2*r2') = for (10) only on the Stiefel manifold, and so level sets of the 
error D{Q) — I are in general not preserved, the error equation being 

D' ^ -2?ft{Dn*An). (22) 

Davey's method, on the other hand, is derived precisely from the global re- 
quirement fl*fl' — 0, and so preserves all level sets of D, with associated error 
equation 

D' = 0. (23) 

That is, (21) corrects for spurious growth modes of (10) in directions normal 
to the Stiefel manifold. 

Remark 6 With the corresponding modification 7' = trace {Q^{A — A±)Q)j 
of (10), (21) gives an alternative implementation of the full polar- coordinate 
method, valid on or off the Stiefel manifold. 

An alternative stabilization of (10) is the damped Drury method 

= (7 - Qn*)A{x, A)n + cfi(/ - , (24) 
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suggested by [11], with c > chosen sufficiently large with respect to the 
matrix norm of A. This evidently agrees with (9)-(10) on the manifold 
Q := {Q : Q*n = I}, but, thanks to the new penalty term on the right-hand 
side of the Q equation, has the additional favorable property that Q is not 
only invariant but attracting under the flow. For, defining D{fl) — — I as 
above, we have the error equation 

D' = -2c(/ + D)D-2 ^{Dn*An) . (25) 

Observing that Q and D have matrix norm one, we obtain for c > \\A\\ that 

D' < -2cD -2^{Dn*An) < 0. 

Remcirk 7 A brief calculation reveals that the stabilizing term fl{I — fl*fl) is 
the steepest descent/ gradient flow for the Frohenius norm squared of D. In this 
sense, it is the optimal, least-squares correction, corresponding approximately 
with orthogonal projection onto the Stiefel manifold D = 0. 

Though (24) in principle remains stable in numerically sensitive regimes where 
(21) may fail, in practice, the large value of c that is required for stability 
introduces numerical stiffness imposing unreasonable restrictions on mesh size; 
see Section 4. A more attractive alternative in this situation is the damped 
Davey method 

= (7 - nn+)An + cn{i - (26) 

obtained by combining the above two stabilization techniques, which is glob- 
ally exponentially attracting for any c > 0, with error equation 

D' ^ -2c{I + D)D. (27) 

This would be our own suggestion in numerically sensitive regimes. 

Alternatively, one might employ the damped Bridges-Reich method 

n' = -2Q[{i - Qn*)A{x, x)]n + cn{i - (28) 

proposed in [7], for which the Stiefel manifold Q likewise is attracting for any 
c > 0. This agrees with (10) on Q, since 

[(/ - nn*)A]*fi = A*{i - nn*)fi = o; 

moreover, the undamped Bridges-Reich method obtained by setting c — has a 
skew-Hermitian form favorable for geometric integrators preserving Lie group 
structure [7] . On the other hand, off of the Stiefel manifold, these methods do 
not preserve the desired subspace, since they are not of the canonical form 

n' ^An + QB, (29) 
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where B — a'a~^, obtained by multiplying (7) on the right by a~^; thus, 
they are in some sense trading errors in the normal direction for new errors 
tangential to the Stiefel manifold. 

Finally, we mention the polar factorization method, a discrete orthogonaliza- 
tion method introduced by Higham [17], based on the observation that the 
factor Q in the polar factorization W = QH of a matrix W, Q orthonormal 
and H symmetric, is the closest point (in usual Frobenius matrix norm) to 
W on the Stiefel manifold. In this split-step method, single steps of an ex- 
plicit scheme approximating the original, linear flow W — AW are alternated 
with projection via polar-factorization onto the Stiefel manifold. Though we 
did not test this scheme, Higham [17] reported superior performance in nu- 
merically sensitive situations, as compared to geometric integrator-based con- 
tinuous orthogonalization methods. For further discussion of the continuous 
orthogonalization method and its variants, see [3,5,7] and references therein. 

Remark 8 In applying Higham 's method to Evans function computations, we 
expect that it is important to substitute for W = AW the asymptotically triv- 
ial rescaled ODE, W' — {A — trace {Q,*AQ,{±oo))I)W , similarly as for the 
exterior-product method [8,9,10] or the radial equation in (12). In this ap- 
proach, 7 is obtained as the product Tlj det Hj of determinants of symmetric 
factors H in the projection steps W = QH Q. Alternatively, one could inte- 
grate Drury's ODE (10) for the ODE step and track by the usual, continuous 
ODE (12). 

Remark 9 An analytic variant of Davey's method is the bi-orthogonal method 

n' = {I -Q{Q*Q)-Hr)AQ, 

n' ^ -{I - n{n*n)-^n*)A*n, ^ ^ 

introduced by Bridges [5]. Here, the matrix fl, and likewise fl* , is analytic along 
with its kxk minors. This scheme has the property of Davey's method that level 
sets ofQ*Q are preserved; in particular, Q*Q = Ikxk is an invariant manifold. 
However, analyticity of fl and fl* is incompatible with fl = fl, and so fl*fl ~ 
I kxk does not enforce good conditioning of fi. Apparently, the requirement 
of analyticity of individual columns is overly rigid for purposes of numerical 
stabilization; the advantage of the polar coordinate method is the flexibility to 
choose optimally the angular evolution equation. 

In our numerical experiments, the best and fastest performance was exhib- 
ited by the original, undamped Drury method (10); see Section 4. Indeed, our 

results indicate that the examples considered lie in the numerically insensi- 
tive regime for which normal instabilities remain small. However, the stability 
issues discussed above may be important for other, more numerically taxing 
problems. 
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Radial equation. Likewise, since they arc scalar quantities, we could more 
stably solve for radial variables 7±(0) by quadrature, using formulae 

trace (Q* Afl(x)~fl* A-fl-)dx ~ / \ 

/-v"y - - -°° ' ' - 7(-oo), ^^^^ 

7+(0) = e-fo^'^"^'' i^*A^(^)-^'+^+^+^d''=^i^+oo). 

However, in practice, this does not seem necessary, as the Q-equation appears 
to be the limiting factor determining accuracy/ allowable mesh size in (12). The 
rescaling (11) on the other hand is critical for good numerical performance, 
as is the analogous rescaling for the exterior-product method [8,9,10], giving a 
speedup on the order of the spectral radius of A{x, A). 



3.2 Continuous initialization 



Though the initializing step at plus and minus spatial infinity is a one-time 
cost, hence essentially negligible, we point out that this step too may be carried 
out more efficiently by an evolution scheme in A parallel to the one carried 
out in X in Section 2.1. Defining W- :— oo,A), 7_ := 7(— oo,A), and 

Q :— [P', P], recall that the A evolution of W- is again given by a finear ODE, 

WL = Q{X)W_. 



Thus, we may apply exactly the same steps as in the previous section to obtain 
a well-conditioned C'' evolution scheme for (n_,7_) of 



n'_ = {i- n_n*_)Q{x)n^ + c{\\Q\\)n_{i - oio. 

7^ = trace (OigO_)7_, 
initializing (7-0, ^^-o) = (1, ^-0) at some base point Aq. 
With this modification, equations (31) simplify to 

~ r" trace (Q*ACl(x)—Q* A-d-)dx ~ / \ 

7_(0) = e-'-'^ ^ V / - 

7+(0) = e/o^°° "^'^'^^ i^*A^(^)-^*+A+n+)dx~^j^^y 



(32) 



(33) 



Note that the above calculation gives the interesting information of the wind- 
ing number of 7_ over one circuit around a closed A-contour, which is not 
necessarily zero, or even an integer. 

Remark 10 A similar, but more complicated scheme could be used to restore 
analyticity in general collocation methods, by further tracking in x the values 
of 7 relating the bases obtained by Kato 's algorithm applied in variable A with 
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respect to orthogonal projection; for, the associated x -evolution depends only on 
the associated subspace, which can be well- approximated, rather than exterior 
product or directions of individual solutions, which cannot. Combined with 
the above computation relating analytic bases at x — oo to those obtained 
through orthogonal projection, we obtain an analytic scheme for which the 
only information required is knowledge (to reasonable tolerance) of subspaces 
at each x. Of course, as pointed out in the introduction, a simpler solution 
would be to use a collocation scheme based on the algorithm of this paper, for 
which no such corrections would be necessary. 



4 Numerical compcirisons 



In this section, we compare our new method for Evans function computa- 
tion with the exterior- product method described in Section 1, and afterward, 
briefly, with various alternative continuous orthogonahzation methods substi- 
tuted in the angular equation. 

As a test system, we consider solitary waves of the "good" Boussinesq equation 

"^tt — "^xx "^xxxx ("^ )xxj (^'^) 

which have the form — x — st) 



(35) 



where the wave speed s satisfies \s\ < 1. We remark that this is the same 
system studied in [2], and is known to be stable when \ < \s\ < 1 and 
unstable when \s\ < 1/2. 



By linearizing (34) about the traveling wave (35), we arrive at the eigenvalue 
problem 

X\ - 2s\u' = (1 - s^)u" - u"" - {2uu)", (36) 
which can be written as a first-order system (2), where 



A{x,X) 









1 








1 





^-A^ - 2u^^ 2Xs - Au^ (1 - s2) _ 2u Oj 



(37) 



For A in the right-hand plane, we have that (37) spectrally separates into two 
growth and two decay modes, that is. A; = 2. We remark that this model is 
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a great test-case for Evans function computation since it captures one of the 
chief numerical obstacles, that is, overcoming multi-mode growth and decay 
in the left and right subspaces, respectively. It is precisely this difficulty that 
motivated the development and use of the exterior-product and continuous 
orthogonalization methods in the first place. 



4.1 Algorithms 



Using the exterior-product method, we lift A{x, A) into exterior-product space 
to get 



A(2)(x,A) 



1 

2As - Au^ (1 - s^) _ 


A2 + 2u^^ 
A2 + 2u^ 




1 

2u 





1 





fl - 






1 
1 

2u 



-2Xs + 4u^ 









1 





(38) 



To maintain analyticity, we use Kato's method (Proposition 3, see also [10,6]) 
for analytically choosing the (simple) eigenvectors r+(— M) and r_(M), which 
correspond to the largest growth and decay modes, respectively, at the nu- 
merical approximates of negative and positive infinity, ±M (we used M = 8), 
where the cigenprojection P is obtained as {l*r)^^rl* for any left and right 
eigenvectors / and r for the eigenvalue of (38) of smallest (resp. largest) real 
part, obtained through a standard matrix routine. To integrate (15) we use 
Euler's first-order method for convenience. We then evolve these vectors from 
X — ±M to a; = using a standard numerical ODE solver (RKF45) and 
compute the Evans function via the wedge product at x = as in (3). 

With our new method, we similarly evolve, analytically in A, the eigenvectors 
at the numerical end states ±M so that we can determine the det(Q*l^='=) 
multipliers in (31). We likewise use Kato's method, where the cigenprojection 
P is obtained as {L* R)^^ RL* for any matrices L and R with columns forming 
orthonormal bases for the left and right stable (resp. unstable) subspace of A, 
obtained by the singular- value decomposition (SVD). For a more efficient, but 
slightly more complicated algorithm, see [18]. 

We likewise determine orthonormal bases r2(±oo) at each A for the stable 
(resp. unstable) subspace of ^4+ (resp. AJ) via the SVD and initialize 7(±M) 
through (20). Finally, we evolve (^^,7) from x = ±M to x = (via RKF45) 
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0.05 0.1 0.15 0.2 -4 -2 2 4 6 8 10 

(a) (b) X 10-^ 

Fig. 1. We evaluate the Evans function of the "good" Boussinesq system for the 
unstable pulse having wave speed s = 0.4. We compute (a) the graph of the Evans 
function along real axis from A = to A = 0.2 and (6) the image of the closed 
contour r{t) = 0.16 + O.OSe^^^**, where 0<t<l. 



and compute the Evans function as the Grammian determinant (5), as de- 
scribed in Section 2. 



4-2 Results 



In Figure 1, we compute the Evans function for the unstable pulse with s = 0.4. 
We perform two Evans function computations for each contour, using our new 
method and the exterior-product method. We then compare results. In Figure 
1(a), we compute the graph of the Evans function along the real axis from 
A = to A = 0.2. From this wc see that the graph crosses through zero 
near A = 0.155, indicating an unstable eigenvalue there. We remark that 
the graphs for both methods were plotted, however since they are virtually 
indistinguishable, it only appears as though one curve is present. In Figure 
1(b), we compute the Evans function about the contour r(t) = 0.16 + 0.05e^'^**, 
where < t < 1, and plot its image. The interior of this contour contains 
the above mentioned unstable eigenvalue, and so we expect the origin to be 
contained in the interior of the image, as it is. This is a second verification of 
the unstable eigenvalue. In this second test, both methods are likewise graphed 
and overlap to the point that they are also indistinguishable. Indeed the two 
graphs differ by an absolute difference of 1.4 x 10^^ and a relative difference of 
4.6 X 10~^. In both old and new shooting algorithms, the absolute and relative 
tolerances are set to 10~^ and 10~^, respectively. 



16 



4-3 Performance 



Function evaluation for (12) requires as many as five matrix-matrix multipli- 
cations, whereas the exterior-product method requires a single matrix-vector 
multiplication. Because of this, our new method is actually slower for the 
"good" Boussinesq, which is a relatively small system. To leading order, the 
operational count for our method grows as 2kn^ + Sk'^n. By contrast, the 

exterior-product method grows as , and is faster than our method when 
n — 4 and k — 2. Indeed for k ~ n/2, we expect the break-even point to be 
at around n — 6. 

We remark that A^''\x,\) becomes sparse as n gets large (k ~ n/2). We 
can show that the number of non-zero entries (and hence operations for a 
sparse matrix- vector evaluation) is exactly {k{n — k) + l)(fc)- Even though 
a sparse-matrix package would offer substantial computational savings over a 
non-sparse one, the size is still prohibitive for large systems. In fact, with a 
sparse improvement, the break-even point for function evaluation relative to 
our method is only extended only to around n = 8. 



44 Other Methods 

We compare the different continuous orthogonalization methods discussed in 
Section 3.1 and examine overall performance and accuracy. Specifically, we 
measure (i) how closely the trajectories "stick" to the Stiefel manifold by 
computing the distance (in Frobcnius norm) of i7 at a; = to the Stiefel 
manifold, {ii) the number of mesh points that are needed for our adaptive 
ODE solver to maintain tolerance (AbsTol=le-8 and RelTol=le-6), {in) the 
run-time needed to compute the Evans function, and finally (iv) the absolute 
and relative errors compared with the exterior-product method taken to high 
accuracy. See Table 1 for the results. 

According to the data, the overall best method for our problem is the origi- 
nal undamped Drury method (12). It is only slightly less accurate than the 
damped Davey method (26) with (c = 5) but is about 25% faster as fewer 
operations are needed. We remark that all three methods become both less 
accurate and (considerably) more time consuming, due to numerical stiffness, 
as the damping constant c gets large. 

We also remark that, while the Bridges-Reich method (28) did not perform 
as well as either the Drury or Davey methods, the undamped method was 
actually designed for use in geometric integrators, and the damping term was 
only introduced to facilitate less sophisticated numerical integration methods. 
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Table 1 



A comparision of methods: The third column "Mesh" correpsonds to the (typical) 
number of mesh points integrating from x = itM to a; = 0. The "Time" column 
measures the run-time for computing the 20 Evans function values along the contour 
7(t) = 0.16 + 40i + .15e2'^**. The last two columns measure the absolute and relative 
difference compared to the values returned using exterior-product method to high 
accuracy. 

Indeed the use of geometric integrators seems like a very good direction to 
explore for Evans function computation, and we intend on exploring this issue 
in future work. 



5 Concluding renictrks 

The numerical results show that the above-described polar-coordinate algo- 
rithm is indeed feasible, and compares favorably in performance even for low- 
dimensional systems to the exterior-product method that is the current stan- 
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dard. Due to better dimensional scaling, it should work well also for large 
systems, whereas the exterior-product method quickly becomes dimensionally 
infeasible. As test problems of medium size, we intend next to investigate 
stability of traveling waves in magnetohydrodynamics and detonation with 
large number of reactant species {n ~ 15). A longer term project might be 
to implement polar-coordinate based collocation methods for extremely large 
systems. 

It seems worth emphasizing a philosophical point associated with the new al- 
gorithm that is simple but possibly of wider use, concerning the apparently 
conflicting goals of numerical well-conditioning vs. maintaining analyticity. 
Namely, in the context of exterior products, an optimal strategy is to pro- 
jectivize, choosing the most numerically convenient basis for the associated 
subspace (Grasmannian), then correct for the resulting loss of analyticity by 
an appropriate scalar factor. Because scalar equations are always numerically 
well-conditioned, the final step costs essentially nothing. That is, we may re- 
cover analytic continuation of subspaces by a simple, well-conditioned post- 
processing step appended to a more standard continuation routine. 

Finally, we point out that our experiments indicate that essentially any exist- 
ing continuous orthogonalization methods should suffice for stability compu- 
tations in low- and medium-frequency regimes relevant for sectorial operators- 
in particular, those arising in the reaction-diffusion context of [23] . For high- 
frequency computations in numerically sensitive situations, it may be neces- 
sary to use the more stable schemes described in Section 3.1. Here, we may 
draw on the resource of the very active community studying numerical con- 
tinuous orthogonalization. 
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